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Self-organization of active matter as well as driven granular matter in non-equilibrium dynamical 
states has attracted considerable attention not only from the fundamental and application view¬ 
points but also as a model to understand the occurrence of such phenomena in nature. These 
systems share common features originating from their intrinsically out-of-equilibrium nature. It 
remains elusive how energy dissipation affects the state selection in such non-equilibrium states. As 
a simple model system, we consider a non-equilibrium stationary state maintained by continuous 
energy input, relevant to industrial processing of granular materials by vibration and/or flow. More 
specifically, we experimentally study roles of dissipation in self-organization of a driven granular 
particle monolayer. We find that the introduction of strong inelasticity entirely changes the nature 
of the liquid-solid transition from two-step (nearly) continuous transitions (liquid-hexatic-solid) to 
a strongly discontinuous first-order-like one (liquid-solid), where the two phases with different effec¬ 
tive temperatures can coexist, unlike thermal systems, under a balance between energy input and 
dissipation. Our finding indicates a pivotal role of energy dissipation and suggests a novel principle 
in the self-organization of systems far from equilibrium. A similar principle may apply to active 
matter, which is another important class of out-of-equilibrium systems. On noting that interaction 
forces in active matter, and particularly in living systems, are often non-conservative and dissipative, 
our finding may also shed new light on the state selection in these systems. 

PACS numbers: 64.70.D-, 45.70.-n, 64.70.ps, 05.70.Ln 


I. INTRODUCTION 

Despite the fact that self-organization of a system in 
an out-of-equilibrium state plays a crucial role in dynam¬ 
ical structural formation in nature, physical principles 
behind such phenomena have remained elusive. Active 
matter [T] and driven granular matter nm are two im¬ 
portant classes of out-of-equilibrium systems. They share 
an intrinsic out-of-equilibrium nature, and the only basic 
difference is that the energy is injected locally for the ac¬ 
tive systems whereas globally for the granular systems [T] . 
This global nature of energy input makes granular mat¬ 
ter physically simpler than active matter. Thus, granular 
matter is not only important for its own sake, but also re¬ 
garded as a model for understanding the physics of active 
matter. 

Granular matter is an important class of materials, dis¬ 
tinct from thermal systems since the thermal energy is 
negligible for its description. Granular matter is ubiqui¬ 
tous in nature and its dynamical self-organization always 
takes place in a strongly non-equilibrium situation as in 
active matter, since energy input is essential for its oc¬ 
currence US]- Its statistical yet athermal nature makes 
the physical description extremely difhcult. From an ex¬ 
perimental point of view, the control of self-organization 
of granular matter is also a difficult task. However, a no¬ 
table exception is a dynamic steady state, maintained by 
the balance between energy input and dissipation, which 
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allows US to perform well-controlled experiments. The 
most idealized system may be a quasi-two-dimensional 
(2D) driven granular particle monolayer, where spherical 
particles are confined between two parallel plates whose 
gap is narrow enough to avoid particle overlap along the 
vertical direction and energy is injected by vertically vi¬ 
brating plates. This system allows us to access all phase- 
space information at the particle level. So the phase be¬ 
havior of such a monolayer particle system has played a 
crucial role in our understanding of the fundamental na¬ 
ture of self-organization in a system far from equilibrium. 

This vibrated monolayer particle system has also at¬ 
tracted considerable attention for its connections with 
fundamental problems in the field of condensed matter 
and statistical physics [5]. The liquid-solid transition in 
a 2D disk system, the thermodynamic counterpart of a 
vibrated monolayer, has been a hot topic since the discov¬ 
ery of the liquid-solid transition for hard disks by Alder 
and Wainwright [^. Two-dimensional particle systems 
cannot crystallize at finite temperature due to significant 
fluctuation effects associated with the low dimensional¬ 
ity, yet the above work shows that they may still form 
solids. There is a long-standing debate [7] on the nature 
of this transition for a system of the simplest interparti¬ 
cle interaction, hard disks. One scenario is that ordering 
takes place via two steps, liquid-to-hexatic and hexatic- 
to-solid transitions, now widely known as the Kosterlitz- 
Thouless-Halperin-Nelson-Young (KTHNY) scenario [51 
IHHID]. Here each transition is continuous. The other 
is that ordering takes place in one step via a first-order 
liquid-solid transition m- There have been hot debates 
on which is the relevant scenario. Very recently, it was 
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shown m that the transition actually takes place by a 
scenario different from both of them: It occurs with two 
steps as the KTHNY scenario suggests, but the first tran¬ 
sition is not continuous but weakly discontinuous. How¬ 
ever, the first-order nature of the liquid-hexatic transi¬ 
tion is very weak, and the transition roughly obeys the 
KTHNY scenario. This basic behavior is common to 
other systems including particles interacting with soft 
repulsive potentials [CTfTH] and those with attractive 
potentials such as the Lennard-Jones potential [17j . al¬ 
though it has recently be shown that the nature of the 
transitions depends on the softness of the potential in 
a delicate manner m- Monolayer granular matter has 
provided a model experimental system to study this fun¬ 
damental problem. 

Some time ago, careful experiments were made on the 
athermal counterpart of the above system. It was shown 
by Shattuck and his coworkers that a driven monolayer 
particle system continuously transforms from a liquid to 
an intermediate hexatic, and then to a solid phase, with 
an increase in the particle area fraction cj) under a con¬ 
stant r m (see Sec. H for the precise definitions of T 
and (j)). A similar meting transition behavior was also 
observed by Olasfen and Urbach when increasing the di¬ 
mensionless acceleration T at a fixed particle area frac¬ 
tion (j) for a granular quasi-monolayer m- However, it 
was shown that the thickness of the cell h, which is 1.6 
times of the particle diameter d, plays a crucial role in 
the transition: height fluctuations of particles may be a 
source of disorder. The increase of their amplitude with 
an increase in the vibration amplitude, or T, increases 
the number density of defects, eventually leading to the 
melting of the solid phase. Thus, the mechanism may 
be essentially different from the former example, which 
does not involve any significant hight fluctuations due to 
a strong 2D confinement. 

The former liquid-solid transition behavior as a func¬ 
tion of (j) [TS] obeys the KTHNY scenario [SJ [5HTU] , al¬ 
though the liquid-to-hexatic transition may be weakly 
first-order da [16]. This study suggests that a quasi- 
2D driven granular system behaves very similarly to its 
thermal counterpart. A similar conclusion was also de¬ 
rived for glass-transition-like phenomena of driven binary 
mixtures [20] and polydisperse systems [211 Ea. The en¬ 
ergy injected by mechanical vibration is converted to the 
kinetic energy of a system and the effective (granular) 
temperature T* is defined by this kinetic energy. If this 
is high enough to overcome the gravity and the energy 
loss originating from the friction and inelastic collisions 
with the container |23j . we may approximately regard 
the system as a thermal equilibrium system as long as 
interparticle collisions are almost elastic. However, the 
exact mechanism responsible for this apparent thermal 
behavior of an athermal system has remained elusive. 
We note that these experiments were performed by using 
rather elastic balls such as steel balls. Then, the natural 
question to ask next is how the nature of the liquid-solid 
transition is affected by the inelasticity of collisions, or 


internal dissipation. 

There were pioneering works on liquid-solid transi¬ 
tions of quasi-2D granular systems . These studies 

showed interesting monolayer liquid-bilayer solid coexis¬ 
tence in a nonequilibrium steady state, in which the two 
phases have different granular temperatures. These are 
intriguing examples of phase ordering, more generally 
self-organization, in a dynamic steady state of a non¬ 
equilibrium open system, maintained by the balance be¬ 
tween energy input and dissipation. 

Here we study the effect of energy dissipation on liquid- 
solid transition of a quasi-2D driven granular matter by 
comparing the behaviors of steel and rubber ball systems. 
We note that steel and rubber balls have differences in 
the restitution coefficient, the friction coefficient, and the 
elastic properties. Among these, the difference in elastic 
properties may be less important compared to the others 
because forces acting upon interparticle collisions are too 
weak to cause non-linear deformation of the balls and it 
is known that the 2D melting behavior is not affected 
by the softness of the interaction potential [16] . A situa¬ 
tion we consider is a single layer of monodisperse spheres, 
which is vibrated between two horizontal plates [30j , un¬ 
like the above-mentioned previous works where bilayer 
formation is allowed [24l - l2^ . The control parameters in 
our experiments are the energy input characterized by 
a dimensionless acceleration T and the area fraction of 
particles 4>. Here we demonstrate that the dissipation 
due to the inelasticity of collisions and the friction can, if 
they are strong enough, completely break the similarity 
between the athermal and the corresponding thermal sys¬ 
tem and fundamentally change the nature of the liquid- 
solid transition in a monolayer from the KTHNY-like 
continuous transitional behavior |18j to a strongly discon¬ 
tinuous transition. We discuss a physical mechanism re¬ 
sponsible for this unconventional self-organization under 
energy dissipation. Our study reveals a novel mechanism 
leading to the coexistence of two phases with different 
granular temperatures, which is essentially different from 
the mechanism previously found [24H29]. We infer that 
a similar mechanism may be relevant to self-organization 
of active matter including living systems. Dissipative in¬ 
teractions, such as inelastic, hydrodynamic, frictional in¬ 
teractions, may play a crucial role in the state selection. 


II. EXPERIMENTAL 
A. Experimental systems 

Our experimental apparatus is analogous to those used 
in m [T^ and the same as that used in [JT] [22]. A 
monolayer of non-cohesive particles was vibrated sinu¬ 
soidally in the vertical direction by an electromagnetic 
shaker (Labworks ET-I39) with frequency / = 50 Hz. We 
changed the dimensionless acceleration T = A{2TTf)'^/g, 
where A is the amplitude of vibration and g is the grav¬ 
itational acceleration, by controlling A. The cell has a 
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circular shape and its inner diameter is L = 102.5 mm. 
The circular annulus made of duralumin is used as the 
hard side wall. When we make an experiment with a 
soft side wall (see Sec. II.D), the inner side wall of the 
annulus is covered by a silicone rubber sheet of 0.2 mm 
thickness. In our experiments, we use two types of spher¬ 
ical particles: steel and fluorine-rubber balls. The par¬ 
ticle diameter is d = 3.0 mm for both stainless steel and 
fluorine-rubber ball. The fluorine-rubber ball has four 
important characteristics: (i) a non-cohesive character, 
(ii) a large stiffness 1 GPa), which allows us to ignore 
nonlinear deformation upon collision, (iii) a low dynamic 
friction coefficient (< 0.4), and (iv) a low restitution coef¬ 
ficient. The friction coefficient of the fluorine rubber may 
be comparable to that of the steel. However, the coeffi¬ 
cient of rolling friction may be smaller for the steel ball 
than for the rubber ball since the former has a smoother 
surface than the latter. 

The restitution coefficient a is difficult to estimate ac¬ 
curately and it is known to depend on the particle veloc¬ 
ity. Here we measured the height ratio before and after 
a collision with a steel wall by dropping a ball vertically. 
The restitution coefficient a against a steel wall was es¬ 
timated as ^ 0.8 and ~ 0.3 with the variance of O.I for 
the steel and rubber ball respectively. We also estimated 
the a for collision between two rubber balls as ~ 0.1. 
The value of a for steel balls coincides well with the lit¬ 
erature value m- By confining these particles between 
two plates separated by an annulus with a thickness h 
of 4.0mm (i.e., h/d = 1.33), we allowed only quasi-2D 
particle motion, i.e, no bilayer formation. We note that 
in previous similar works on the liquid-solid transition 
of driven granular matter [26l [28] h/d = 1.74 ~ 1.95, 
which allows particles to form bilayers (see Appendix A- 
1). This difference in h/d leads to a crucial difference in 
the nature of a liquid-solid transition between our system 
and theirs, as will be discussed later. 

The top plate is a transparent glass plate so that we 
can observe the particles through it. The bottom dura¬ 
lumin plate is covered with a sandpaper by which the 
vertically oscillating particles are scattered [TH] . The use 
of a wall with surface roughness is expected to randomize 
the horizontal motion and realize Brownian-like motion 
of particles. It was shown [32| that a large enough T leads 
to the velocity distribution function of a nearly Gaussian 
shape for a driven monolayer system. This is because 
the randomization of the energy injection due to particle 
collisions and the wall roughness leads to realization of a 
self-generated effective white bath in a long timescale. 

For steel balls, we obtained the particle coordinates by 
tracking the bright spots at the top of particles reflecting 
illuminating light. For rubber balls, on the other hand, 
we obtain the particle coordinates simply by binarisation 
of particle images. We used a high speed camera (VCC- 
H500, DigiMo Co. Ltd.). The typical frame rate used 
was 100 frame/s and the image resolution was 640x480. 
Occasionally, we used the rate of 500 frame/s with a res¬ 
olution of 640x90. A pixel size corresponds to 0.25 mm. 


The position of each particle was tracked by a particle 
tracking software, which fits a Gaussian function to an 
image of each particle. The detection error is less than 
0.05 mm for a steel ball and 0.1 mm for a rubber ball. 
We measured a structure after attaining a steady state 
(typically after 10 min from the initiation of vibrational 
driving). 

In our system, the transition between a solid-like and 
liquid-like state occurs as a function of the area frac¬ 
tion (f) = {Nd'^)/L'^, where d is the diameter of particle, 
L = 102.5 mm is the inner diameter of the annulus, N is 
the total number of particles in the system. We note that 
this definition of (/ is often used in describing the phase 
transition of (quasi-)2D hard particle systems (see, e.g., 
[18)1. We note that the random close packing and the 
closest packing in 2D occur at about 0.83 and 0.906 re¬ 
spectively. 

B. Finite-size effects 

Because of the rather small size of our system, our 
measurements may suffer from finite-size effects. The 
effects may not be so serious for a discontinuous liquid- 
solid transition in a rubber ball system, since there is no 
diverging lengthscale. On the other hand, the results of a 
steel ball system may suffer from finite-size effects. How¬ 
ever, the maximum correlation length we could attain 
was 9 times the particle diameter d, which is still much 
shorter than the cell diameter (^ Sid). So the finite 
size effects may not be so severe for our results, although 
there might be slight shifts in the transition area frac¬ 
tions. This conclusion is supported by the rather good 
agreement of the measured ^-dependences of the correla¬ 
tion length of hexatic order and translational order ^ 
with the theoretical predictions for an infinite system (see 
below). However, the correlation lengths, and may 
be bound by the cell size very near the hexatic ordering 
point (ph and the solidification point (ps, respectively, and 
the weak discontinuous nature of a liquid-hexatic tran¬ 
sition may be smeared out by the finite-size effects [12) . 
In principle, we can use a cell with a larger diameter, 
but what is most important in our measurements is the 
homogeneity of the vertical vibration amplitude. So we 
used only the rather small cell. We are planning to make 
a large cell with high mechanical rigidity, but we leave 
this for future investigation. 

C. Characterization of structures 

The 2D radial distribution function g{r) was calculated 
as 

which is the ratio of the ensemble average of the number 
density of particles existing in the region r ~ r -|- Ar 
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to the average number density p = N/Li^. Here N is 
the number of particles in the simulation box, whose side 
length is L, and Ar is the increment of r. The spatial 
correlation of translational order is characterized by the 
translational correlation length ^ ioi (j) < (ps as follows: 

En+{g{r) - 1) ^ exp(-r/^). 

where Eri^{f{x)) expresses an operation to extract the 
envelop function of the positive part of a function f(x). 

The hexatic order parameter is measured by i/'e = 
(l/rij) exp (idOjk) for each particle j, where the sum 
runs over the rij nearest-neighbors of particle j, and is 
the angle between the bond — rj and a fixed arbitrary 
axis. Here the nearest-neighbor particles are detected 
by a criterion r < 1.4c?. We confirmed that this crite¬ 
rion provides the same nearest-neighbor identification as 
a method based on Voronoi tessellation. 

Then, the spatial correlation of ip^ is calculated as [5] 

= 2^rArW(JV - 1) S 

The spatial correlation of the bond-orientational order 
can then be characterized by gQ^{r)/g{r). Here the di¬ 
vision by g(r) is to remove the effect of translational or¬ 
dering. 

To characterize the fluctuations of tpQ, we estimate the 
spatial correlation length of ipe, and the susceptibility, 
Xe = (('06 — (06))^)- We obtain the spatial correlation 
length ^6 for (p < (ph from the following relation: 

En^{9l^ir)/g{r)) ^ exp(-r/^6)- 

III. RESULTS 

A. Ordering in a steel ball system npon 
densiflcation. 

Now we show experimental data which provide infor¬ 
mation on the nature of the phase transitions. For driven 
steel balls, it was previously shown that the system trans¬ 
forms from the liquid to solid via the intermediate hex¬ 
atic phase m- For comparing the behavior of steel and 
rubber ball systems on the same ground we performed 
experiments for the two systems under the same experi¬ 
mental conditions. 

Before showing results, we briefly review how the spa¬ 
tial correlations of the positional and bond orientational 
order parameter are predicted to grow with an increase 
in the area fraction 0 in the framework of the KTHNY 
theory [S]. The predicted behaviors are summarized in 
Fig.E The hexatic ordering point <ph and the crystalliza¬ 
tion point (ps are characterized by the power-law decays 
of the bond orientational and positional order parame¬ 
ter, respectively. The former should obey at (ph, 

whereas the latter at (ps- 

We first describe the results of a steel ball system. We 
show the 0-dependence of g{r) and gQ{r)/g{r) for a steel 
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FIG. 1. Phase behavior of 2D hard disks predicted by the 
KTHNY theory. In this table, we summarize how the spatial 
correlations of the positional and bond orientational order 
parameter should increase with an increase in 0. The other 
characteristics of each phase are also shown. 

ball system in Fig. The spatial decays of both quan¬ 
tities change from an exponential to a power-law decay, 
we identify the ordering points by whether the decay is 
slower than the predicted power law decay or not. Be¬ 
cause of the limitation of our system size, the firm con¬ 
firmation of the exponent of asymptotic power-law decay 
of these correlation functions are difficult. But the re¬ 
sults are at least consistent with the prediction of the 
KTHNY scenario (see also below). As shown in Fig. 

En'^{g1^{r)/g{r)) should decay slower than for 

'P ^ 4'h, whereas E'nP' [gfg^ should decay slower than 
for 0 > 0s. In Figs, ^a) and (b), these lines with 
a slope of —1/4 and —1/3 are drawn as guides to estimate 
ph and 0s respectively. 

In this way we determine ph ~ 0.72 and ps ~ 0.735. 
The determinations of ph and ps are further supported 
by the divergence of the correlation length and the sus¬ 
ceptibility of the hexatic order parameter for the former 
and by the divergence of the translational order corre¬ 
lation length for the latter. The slight differences of the 
transition area fractions from those of a hard disk system 
[12] may be due to the quasi-2D nature of a system and 
finite size effects. 

As examples of our analysis, here we show details of 
our fittings of ge{r)/g{r) for the steel ball system, respec¬ 
tively, in Fig. ^ The fittings to g{r) — 1 are basically the 
same. From these analyses, we obtain the 0-dependence 
of the spatial correlation length of the hexatic order 06, 
^6, and that of the translational order, although there 
are considerable errors due to the small system size. The 
growth of these lengths with an increase in p are used to 
confirm the prediction of the KTHNY theory in a quan¬ 
titative manner (see Fig. B^) and (b)). We can also see 
indications that the decays of the envelops of ge{r)/g{r) 
and g{r) — 1 change from exponential to power law, re¬ 
spectively, around ph = 0.72 and ps = 0.735. 

The phase boundaries between liquid, hexatic, and 
solid phases are determined by standard methods used 
to study a thermal system, such as the nature of the de- 
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FIG. 2. 0-dependence of g{r) and ge{r)/g{r) for a steel ball 
system for F = 3.3. To visually confirm a decay slower than 
the asymptotic decay for ordered phases, we plot the data in 
the log-log plots, where only the points having positive val¬ 
ues are displayed, (a) 0-dependence of g{r). The signals are 
0 =0.70 (dark blue), 0.71 (yellow), 0.72 (green), 0.73 (orange), 
and 0.74 (red) from the bottom to the top. The straight line 
has a slope of -1/3. (b) 0-dependence of g%(r)/g{r). The sig¬ 
nals are 0 =0.65 (right blue), 0.70 (dark blue), 0.71 (yellow), 
0.72 (green), 0.73 (orange), and 0.74 (red) from the bottom 
to the top. The straight line has a slope of -1/4. From these 
results, we identify 0 =0.65, 0.70 and 0.71 as a liquid state, 
0 =0.72 and 0.73 as a hexatic state, and 0 =0.74 as a solid 
state. 

cay of the spatial correlation of bond orientational and 
translational order (see Fig. [^. We measured the cor¬ 
relation length of the hexatic order parameter ^g, 
as a function of the area fraction 0. Firstly, the de¬ 
cay of the correlation in the hexatic order parameter 
changes from an exponential to a power-law decay at 
the liquid-hexatic transition point, 4>h (— 0.72); secondly, 
the density correlation changes from an exponential to 
a power-law decay at the hexatic-solid transition point 
0s (= 0.735) (see above). As shown in Fig. Ha), we 
observe the steep divergence of ^g towards 4>h, which 
is consistent with the prediction of the KTHNY theory, 
^6 = ■ We also fit the re¬ 

lation of ^ = ^ 0 exp(c|(l/0s) — (1/0)1“°-^^) to the data 
of although the accuracy of the data is not so high. 
We also observe a sharp peak in the susceptibility of 
X6 = ((V'6 - (V'e))^) at (Fig. He)). In this way, we 
confirm the presence of a liquid, hexatic, and solid state 
and determine the phase-transition compositions (j)h and 
0s. However, the estimation of 0s may not be so accu¬ 
rate since the range of the analysis of g{r) — 1 is limited 
by finite size effects (see Sec. II B). Furthermore, the 
short-range order developing for r < 4d do not allow us 
to access an asymptotic power law decay of g{r) — 1 in a 
wide distance range. To avoid finite size effects, we need 
to perform experiments in a large cell. However, we can 
at least see the very slow algebraic decay of g{r) — l above 

0s ■ 

We note that, for steel balls, we do not observe any 
indication of liquid-solid coexistence, also suggesting the 
(nearly) continuous nature of the two transitions. The 
basic phase ordering behavior of a steel ball system is 
thus fully consistent with the KTHNY scenario, as re- 



FIG. 3. Examples of the analysis of the spatial decay of 
g&{r)/g{r) for a steel ball system for F = 3.3. The blue 
points are used for the fittings. The blue curves are expo¬ 
nential fittings, whereas the orange lines have a slope of -1/4, 
which is expected at 0;,. From the above, we can judge that 
the liquid-to-hexatic transition takes place between 0 = 0.715 
and 0.72. The accuracy of this determination is limited by the 
discreteness of the area fractions we employed. 

ported previously by Shuttack and his coworkers [18]. 


B. Ordering in a rubber ball system upon 
densification. 

Next, we focus on the liquid-solid transition in a rub¬ 
ber ball system. Unlike in the above steel ball system, 
we observe in this system the two-phase coexistence of 
the liquid and solid phases separated by a rather sharp 
interface at a certain range of 0 (see Fig. H. The two- 
phase coexistence can also be clearly seen from the bi- 
modal shape of the probability distribution function of 
the hexatic order parameter ^(^jg), as shown in Fig. 
Below 0i, the spatial correlations of hexatic and trans¬ 
lational order are both short-range and decayed nearly 
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(j) (j) 



<t> ^ 


FIG. 4. Continuous and discontinuous nature of the phase 
transition respectively in steel and rubber ball systems. Blue 
curves and symbols are for a steel ball system, whereas red 
ones are for a rubber ball system, (a) The divergence of the 
correlation length of the hexatic order, (filled circle) to¬ 
wards the hexatic ordering point (ph and that of the transla¬ 
tional order, (open circle) towards the solidification point 
(j >3 for a steel ball system. The solid and dashed curves are 
the prediction of the KTHNY theory (see Fig. Q. (b) The (j)- 
dependence of (filled circle) and ^ (open circle) for a homo¬ 
geneous liquid state of a rubber ball system, (c) The suscep¬ 
tibility of the hexatic order parameter xe for steel (blue filled 
circle) and rubber ball systems (red filled circle). The steep 
increase near is observed for a steel ball system, reflecting 
the nearly continuous nature of the transition, whereas there 
is no such behavior for a rubber ball system, reflecting the 
strong discontinuous nature of the transition, (d) The area 
fraction of the solid phase to the total area. As/A, as a func¬ 
tion of ij>. It starts to increase from <j)L until (j>s in proportional 
to (f> — (j)L, indicating the validity of the lever rule. We also 
confirmed the similar relation between the number of parti¬ 
cles of the solid phase Ns and the total number of particles 
N. This means that each coexistence phase retains the same 
number densities irrespective of tj) in the coexistence region. 
The error bars in (a)-(c) represent the standard deviations of 
the fittings. For all the data, F = 3.3 

exponentially, as shown in Fig. As can be seen in Fig. 
I^b), neither the translational nor orientational correla¬ 
tion length, ^ and exhibit any growth with an increase 
in (j), unlike the case of a steel ball system. There is also 
no increase in the susceptibility around the phase tran¬ 
sition points, as shown in Fig. |^c). Since the analysis 
of the decay of the spatial correlation functions is not 
useful in the coexistence region, we distinguish a liquid 
state, a solid state, and a coexistence state on the basis of 
measurements of the local orientational order, the local 
density (or, area fraction), and the local mobility, which 
has a link to the granular temperature, (see Fig. [^. We 
can see these three different quantities including both 


static and dynamical ones can identify the two-phase co¬ 
existence in a consistent manner, although the analysis 
of an instantaneous structure leads to some inaccuracy 
due to short-time fluctuations of the interface (see below 
on a possible origin). We can also clearly see that the 
fraction of the solid phase monotonically increases with 
an increase in </>■ The area fraction of the solid phase 
determined from the hexatic order parameter is shown 
in Fig. I^d), which indicates the validity of the lever 
rule (see below). A discontinuous first-order-like phase 
transition often accompanies hysteresis and metastabil¬ 
ity. We made all the observation after attaining a steady 
state (typically after 10 min from the initiation of vibra¬ 
tional driving) and the behavior was very reproducible. 
We did not observe any indication of hysteresis mainly 
because measurements are always done after a steady 
state is reached and there is no way to change the area 
fraction continuously. So we identify 4>l and (ps as the 
lower and upper boundary of the solid and liquid phase, 
respectively. 

Here we show examples of of our fittings of g{r) and 
ge{i')/g{r) for the rubber ball system at </> = 0.69 in Fig. 
1^ We find no systematic ^-dependences for g{r) and 
gQ(r)/g{r) (see Fig. Qb)). We can see in Fig. [^that even 
at 0 = 0.69, which is very near (j)]^ = 0.695, the decays of 
density correlation and hexatic order correlation are both 
exponential and much faster than the power law decays 
of exponent —1/3 and —1/4, respectively. 

Above i/s, we see a homogeneous solid phase. For (/l < 
(j) < (j)Si we observe the coexistence of the liquid and the 
solid phases. We confirm that the former has the upper 
bound (j) of the liquid phase, 4>l, and the latter has the 
lower bound (j) of the solid phase, (j)s- We determine 
(j)L = 0.695 and cj)s = 0.775. The area fraction of each 
phase obeys the lever rule in the coexistence region, as 
shown in Fig. Bd). 

The strong discontinuous nature of the transition is 
also confirmed by that fact that there is neither the di¬ 
vergence of ^6 (Fig. ^b)) nor the sharp increase in xe 
(Fig. B c)) at the phase boundary. We also show the 
distribution function of tpQ, in Fig. [^e) and (f). 

We can see a clear bimodal distribution for the coexis¬ 
tence region of (f>, which is another clear indication of the 
liquid-solid coexistence. Such a signature of the strongly 
first-order transition is absent for a system of steel balls. 

We find that the coexistence can be seen not only by 
the local packing symmetry characterised by the hexatic 
order parameter (Fig. Ba)), but also by the area fraction 
(j) (Fig. Bb)) and the in-plane displacement amplitude 
(Fig. Bw). The last point indicates that the effective 
(granular) temperature T* defined by the kinetic energy 
is spatially inhomogeneous for the liquid-solid coexisting 
state. Thus, the solidity is linked not only to the static 
quantities such as the area fraction and the bond orien¬ 
tational and translational order, but also to the effective 
temperature. This is a very unique feature of this non¬ 
equilibrium steady state, which is maintained by con¬ 
tinuous vibrational energy input and dissipation due to 
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FIG. 5. Liquid-solid coexistence for rubber ball systems 
observed at 0=0.70, 0.72, and 0.74, which are between 
0n(=O.695) and 0s(=O.775), as F = 3.3. The two-phase co¬ 
existence can be seen by binarisation by using the following 
three quantities: (a) the hexatic order parameter 06 of each 
particle (here the threshold is chosen as 06*' = 0.6, and the 
pattern is insensitive to the choice for 0 = 0.6 — 0.7 (see Fig. 
§)• Yellow (solid) and black (liquid) particles correspond to 
particles having 06 > 06*' and 06 < 06*', respectively, (b) 
the local density p, which is calculated from the Voronoi area 
of each particle (the threshold 1/p**' = 9.4 mm^). Blue-green 
(solid) and pink (liquid) particles correspond to particles hav¬ 
ing p > p**' and p < p**', respectively, (c) the displacement 
over 10 s of each particle (the threshold (5**'=3.4 mm). Blue 
(solid) and green (liquid) particles correspond to particles 
having <5 < and S > <5**', respectively. We can see that 
the solid regions identified by these three quantities are well 
correlated with each other. We can also see the monotonic 
increase of the solid fraction with an increase in 0. 

inelastic interparticle collisions. This phenomenon has 
some similarity to the inhomogeneization of a granular 
gas due to inelastic collisions m 1 [31 [B]- Coexistence 
of phases with different granular temperatures was also 
reported for a quasi-2D system where bilayer formation 
is allowed [2m29| . however, the underlying mechanism 
may be quite different, as will be discussed later. The 
link between the solidity and the effective temperature 
indicates that the interfacial profile may be related to 
the spatial gradient of not only the area fraction and 
the hexatic order but also the temperature (i.e., the ki¬ 
netic energy). This possible dependence of the interfa¬ 
cial profile on the mobility, or the effective temperature, 
is unique to non-equilibrium open systems and absent in 
thermal equilibrium systems (see, e.g., [33]) • 



% % 



FIG. 6. 0-dependence of the probability distribution func¬ 
tion of the hexatic order parameter 06, P{'ipe), at F = 3.3. 
Blue curves are for a steel ball system, whereas red ones for a 
rubber ball system. For a steel ball system, P(0) always has 
a unimodal shape, whereas for a rubber ball system it has a 
clear bimodal shape for 0 between 0n and 0s. The long tails 
of P(06) toward low 06 in homogeneous solid phases come 
from defects. 



FIG. 7. 0-dependence of g{r) and g%(r)/g{r) for a rubber ball 
system at F = 3.3. (a) 0-dependence of g{r). The signals are 
0 =0.65 (bule), 0.69 (yellow), 0.70 (green), 0.74 (orange), and 
0.80 (red) from the bottom to the top. (b) 0-dependence of 
ge.{r)/g{r). The signals are 0 = 0.65 (bule), 0.69 (yellow), 
0.70 (green), 0.74 (orange), and 0.80 (red) from the bottom 
to the top. 

C. State diagrams. 


On the basis of these results, we draw the phase (more 
strictly, state) diagrams for steel and rubber balls, which 
are shown in Fig. Here we include the dependence of 
the phase behavior on F. In general, we need F larger 
than a critical value F,. to maintain a dynamical steady 
state. Below Fc, the energy injection becomes inhomo¬ 
geneous. The gray region labeled “inhomogeneous exci¬ 
tation” in Fig. ib) is in such an inhomogeneous state. 
This is because a well-defined dynamical steady state can 
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FIG. 8. Examples of the analysis of the spatial decay of g{r) 
and g6{r)/g{r) for a rubber ball system at F = 3.3. Here 
we show the analysis of g{r) (a) and ge(T)/g(r) (b) at (^ = 
0.69, which is very close to (j>L = 0.695. The blue points are 
used for the httings. The blue curves are exponential fittings, 
whereas the red and orange lines have a slope of —1/3 and 
— 1/4 respectively. In a one phase region, the correlation of 
the bond orientational order decays exponentially and there 
is no indication of a power-law decay even near </_l. The 
extracted correlation lengths in this manner are plotted as a 
function of (j) in Fig. Bb). We note that in the two-phase 
coexistence region, this type of analysis is not meaningful. 

be maintained only when the energy injected by vertical 
vibration overcomes the effects of the gravity and the 
energy loss due to inelasticity of collisions of balls and 
the walls [531 IMl ISZ] (see also Appendix A-2 and Ap¬ 
pendix B). We find the critical value of F to maintain a 
dynamical steady state is higher for the rubber ball sys¬ 
tem (Fc ~ 2.7) than for the steel ball system (F,, ^ 1.3), 
reflecting the stronger inelastic nature of collisions of the 
former with the confining plates. We note that the tran¬ 
sition area fractions, (j)h, 4's, (f’L, and 4>s, are all indepen¬ 
dent of F within the accuracy of our measurements. 

The insets of Fig. [^a) and (b) show the local orien¬ 
tational order in the solid phases formed in a steel and 
rubber ball system at </> = 0.78, respectively. We can 
clearly notice that the amount of defects is much larger 
for a rubber ball system than for a steel one. This can 
also be seen in the shape of P{ip6), which has a larger tail 
toward low ijjQ for a rubber ball system than for a steel 
one. This may be just a consequence of the difference in 
the lower stability limit (/> of the solid phase between a 
steel and a rubber ball system, but there might be other 
fundamental reasons. 


D. Dissipation-induced wetting. 

We also note that the solid phase is always formed 
in the middle part of the container far from the side 
wall (see Fig. and Fig. l^a)). This may be ex¬ 
plained by a larger restitution coefficient of particle-side 
wall collision than particle-particle one. The hard side 
wall prefers the liquid phase with high T* . We con¬ 
firm that a softer side wall (covered by a silicone rubber 
film) is statistically more wettable to the solid phase (Fig. 
Mb)), although wetting is rather modest likely because 
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FIG. 9. Phase diagram of driven monolayer sphere systems. 
Filled circles represent state points where the measurements 
were made (blue: disordered liquid; orange: hexatic phase; 
pink: solid; green: a coexistence of liquid and solid; gray: 
no well-defined stationary state due to inhomogeneous exci¬ 
tation). (a) Steel ball systems, which obey the KTHNY-like 
scenario. In this case, the state behavior does not depend 
on F in the range studied, (b) Rubber ball systems, which 
show the distinct discontinuous transition between liquid and 
solid states. In this case, the dynamical steady state is real¬ 
ized only above F = 2 . 7 . This higher threshold value of F is 
presumably due to the inelastic nature of collisions of rubber 
balls with the confining plates. The images in (a) and (b) are 
configurations in the solid phase at (/ = 0.78 for steel and rub¬ 
ber ball systems, respectively. Green, blue, and red particles 
have local hexagonal, pentagonal, and heptagonal structures, 
respectively. The inner part surrounded by the red dashed 
circle has less order (i.e., more defects) for the rubber ball 
system than for the steel one. 

the curved wall inevitably induces elastic distortion of 
the solid phase. Thus, this wetting phenomenon may be 
regarded as dissipation-induced wetting. 


E. Inelasticity-induced demixing. 

Finally, we briefly mention a related interesting phe¬ 
nomenon we observe in a mixture of steel and rubber 
balls. There are some studies on inelasticity-induced 
demixing [SIlEHHIl]; however, there is no example as- 
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FIG. 10. Wall effects on the liquid-solid coexistence, (a) A 
coexistence of the liquid and solid phase confined by a circular 
steel wall for (j) = 0.74. The liquid phase preferentially wets 
the wall, (b) The same as (a) for a softer wall whose surface 
is covered by a silicone rubber film. In this case, the solid 
phase partially wets the wall (see the regions surrounded by 
the red lines). 

sociated with a liquid-solid transition. Here we report 
inelasticity-induced demixing, which is linked to the for¬ 
mation of an ordered solid phase. When we mix a small 
amount of steel balls with rubber balls, we observe steel 
balls are completely expelled from the solid region and in¬ 
cluded in the liquid region (see Fig. [TT|(a)). We note that 
the solid phase is made of only rubber balls: dissipation- 
induced demixing. Here we stress that in such a mixture 
the energy input rate to steel balls is higher than that 
to rubber balls. On noting that steel balls have a larger 
kinetic energy due to its larger restitution coefficient, the 
above argument naturally explains the preferential inclu¬ 
sion of steel balls into the liquid phase. If the fraction 
of steel balls is too large, the coexistence conditions can¬ 
not be satisfied simultaneously and thus the formation of 
the macroscopic ordered solid phase is largely prohibited 
even at the same area fraction 4>- As a result, a rather 
random mixed liquid state is realized (see Fig. ©b)). 

IV. DISCUSSION 

A. Comparison of our study with previous works 
on driven quasi-2D systems. 

Here we consider the relationship of our results with 
those of pioneering works on phase orderings in quasi- 
2D driven granular systems , in which bilayer for¬ 

mation is allowed. We note that there is a crucial dif¬ 
ference in the dependence of the phase behavior on the 
strength of inelasticity: For systems, where a bilayer can 
be formed, ordering (towards bilayer solid) is suppressed 
by inelasticity [33 [33 [33]. We emphasize that this is 
completely the opposite to our case: inelasticity helps 
ordering (towards monolayer solid). See Appendix A-1 
for a more detailed discussion including an intuitive ex¬ 
planation on the cause of the difference. We also note 
that the monolayer-bilayer formation can take place as 
a function of F [33 03] : whereas our transition cannot 
be induced by changing F as we can see from the phase 


FIG. 11. Phase behavior of mixtures of rubber and steel balls, 
(a) Snapshot of a 1:7 mixture of steel and rubber balls (the 
total area fraction of <() = 0.72) driven at F = 3.3. Here we can 
see a coexistence of the solid phase purely made of rubber balls 
and the liquid phase which is a mixture of steel and rubber 
balls, (b) Snapshot of a 1:1 mixture of steel and rubber balls 
(the total area fraction of ()> = 0.72) driven at F = 3.3. The 
color of the outer shell of each particle represents the type of 
the particle (blue: steel ball; red: rubber ball), whereas that 
of the inner core the degree of the hexatic order (black: low 
i/>6; yellow: high 

diagrams in Fig. Near Fc, we see a crossover from 
inhomogeneous excitation to homogeneous one, but this 
transition as a function of F has an origin essentially 
different from the transition as a function of (p (see Ap¬ 
pendix A-2 and B for the details). The above-mentioned 
fundamental differences between our work and the previ¬ 
ous studies strongly indicate that the underlying physics 
is essentially different between the two cases even on a 
qualitative level. In Appendix A-3, we also mention an¬ 
other type of inhomogeneization in driven inelastic par¬ 
ticles, gas-liquid coexistence [44] . 

Finally, we mention the 2D melting of a solid phase 
with an increase in F m- As described in the introduc¬ 
tion, this melting is induced by the increase in the defect 
density, which is caused by the increase in hight fluctu¬ 
ations of particles. Since our 2D constraint /i/d ~ 1.3 is 
much stronger than the one used in this work (h/d = 1.6), 
we did not see any indication of F-induced melting. It is 
interesting to study in the future how and at which h/d 
the behavior changes. 

B. Roles of inelastisity and friction in 
phase-ordering behavior. 

First we qualitatively consider what is the origin of 
the difference in the phase-transition behavior between 
the steel and rubber ball systems. Steel and rubber balls 
differ not only in the inelasticity and friction upon in¬ 
terparticle collisions but also in the softness of interpar¬ 
ticle interactions. It is rather well established that the 
KTHNY-like behavior, or a transition from liquid to solid 
via the intermediate hexatic phase, is very robust for 2D 
‘thermal’ systems, irrespective of the types of interpar¬ 
ticle interactions: Essentially the same behavior was ob¬ 
served not only for hard disks but also for particles inter¬ 
acting with soft repulsive potentials and attractive (e.g.. 
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Lennard-Jones) potentials (see, e.g., Refs. [TMT^ '). al¬ 
though it has recently been shown that the liquid-hexatic 
transition depends delicately on the softness of the inter¬ 
action [TB]. This robustness of the KTHNY-like behav¬ 
ior irrespective of the nature of interactions is a natu¬ 
ral consequence of the fact that hexatic order is formed 
primarily by geometrical packing effects: When hard- 
sphere-like particles are packed, this local configuration 
is entropically favored. This suggests that the nature of 
the liquid-solid transition between the steel and rubber 
ball systems is not controlled by the softness of particles, 
but by the dissipation due to the inelasticity of particles 
and the friction. 

Now we consider how inelasticity and friction can 
change the nature of the liquid-solid transition. Steel 
balls are not perfectly elastic but the system still be¬ 
haves like a thermal system. On the other hand, the 
behavior of the rubber ball system is distinctly differ¬ 
ent from the behavior of its thermal counterpart. So the 
question can be rephrased as what physical mechanism 
controls the transition from apparently thermal to ather- 
mal behavior when we change the degree of inelasticity 
and friction from steel to rubber balls. For a vertically vi¬ 
brated monolayer system, the particle-(top and bottom) 
wall collision frequency, fpw, is higher than the particle- 
particle collision frequency, fpp, as long as T > T^. We 
estimated fpw as roughly 100 Hz by illuminating light 
from a low angle from the horizontal plane and fpp as 
roughly 20-50 Hz, depending upon on the local (j), by us¬ 
ing normal illumination from the top, with a fast camera. 
For an inelastic particle system to apparently behave like 
an elastic thermal system, the energy dissipated by each 
particle-particle collision should be fully recovered by en¬ 
ergy input through particle-wall collisions before the next 
collision with a particle takes place. As long as this con¬ 
dition is satisfied, even a dissipative system apparently 
behaves like a thermal system, as our steel ball system 
does. 

The athermal nature is accompanied by inhomo- 
geneization of a system and the resulting coexistence of 
two phases with different effective granular temperatures. 
Such two-phase coexistence characteristic to athermal 
systems was already discovered for a quasi-2D driven 
granular system, in which bilayer formation is allowed un¬ 
like our strictly monolayer system [MH^ . In such a case, 
the different granular temperatures can be explained by 
the difference between monolayer excitation and bilayer 
excitation (e.g., differences in effective mass and type of 
excitation). However, this mechanism cannot explain the 
two-phase coexistence we found in a driven monolayer 
system. Below we consider the origin of inhomogeneiza- 
tion and the resulting coexistence of two phases with dif¬ 
ferent granular temperatures on a qualitative level. 

Here we consider a gedanken experiment, where the 
inelasticity is gradually increased in a continuous man¬ 
ner. With an increase in the inelasticity of particles, 
the perfect recovery before the next interparticle colli¬ 
sion becomes more and more difficult. With an increase 


in the particle density, or 0, the inter-particle collision 
frequency fpp, which controls the rate of energy dissipa¬ 
tion, increases but with a rather constant particle-wall 
collision frequency fpw, which makes this recovery pro¬ 
cess less efficient. The energy dissipation is due to both 
the inelastic and the frictional nature of interparticle col¬ 
lisions. We believe that this crossover from the per¬ 
fect to imperfect recovery of the kinetic energy during 
1/fpp with an increase in the particle density destabi¬ 
lizes a non-equilibrium steady state with spatially ho¬ 
mogeneous density and granular temperature, leading to 
the transition from an apparently thermal to strongly 
athermal behavior. Now we qualitatively consider how 
the imperfect recovery makes a system inhomogeneous: 
In such a situation, the degree of recovery of the kinetic 
energy of a particle should decrease with an increase in 
local (j) since fpp increases with an increase in (jj. Thus, 
higher density regions formed by spontaneous fluctua¬ 
tions should transiently have lower kinetic energy than 
lower density regions, which leads to lower pressure in 
the former. This mechanism is essentially the same as 
that of clustering instability in dissipative gases [331 • 
the higher density regions are further compressed, re¬ 
sulting in the enhancement of density fluctuations. How¬ 
ever, an increase in density eventually causes the increase 
in the inter-particle collision frequency and thus the in¬ 
crease in pressure. When the horizontal pressure is spa¬ 
tially homogeneized, this development of density fluctu¬ 
ations stops, ending in two-phase coexistence. In this 
way the system finally reaches a steady state, where the 
energy input rate and the dissipation rate is balanced 
while keeping two phases with different (j) and granular 
temperature T*. Thus, a high enough density region and 
a low enough density region should coexist in a steady 
state for a certain range of the average (f>. The aver¬ 
age horizontal pressure in the two phases should be the 
same since the mechanical force balance is to be satisfied 
across the domain interface. Because of its mechanical 
nature, this equal pressure condition may be hold even 
in an athermal condition (see below, however, on the ef¬ 
fect of interparticle friction). However, the fluctuations 
of pressure may be much larger than in a thermal sys¬ 
tem, since the particle velocity is determined not only 
by interparticle collisions but also by particle-wall colli¬ 
sions. The another constraint comes from the condition 
for a steady state: the balance between the energy input 
rate and the dissipation rate should be balanced in each 
phase. The low density liquid phase is characterized by 
low fpp, which means low energy dissipation and results 
in high T*, whereas the high density solid phase by high 
fpp, which results in low T*. This situation is realized 
under a homogeneous energy input, or spatially homoge¬ 
neous fpui- On the basis of this physical picture, below 
we consider a principle behind the two-phase coexistence 
on a phenomenological level. 

Here we note that it is interesting to study the de¬ 
pendence of the phase-transition behavior as a function 
of inelasticity. However, a precise control of particle in- 
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elasticity is not easy in experiments, and thus numerical 
simulations may be more promising. 

C. How to approach the problem 

Before discussing a principle behind the two-phase co¬ 
existing more theoretically, here we consider which theo¬ 
retical framework is suitable for the description of what 
we discussed above. One candidate is a hydrodynamic 
theory for a confined granular system, which was recently 
developed by Brito et al. [IS]. It was shown that under 
a situation where the thermostat does not inject momen¬ 
tum but only energy, the equations for the conserved den¬ 
sity field and momentum density are the continuum equa¬ 
tion and the Navier-Stokes one, as in the case of usual 
fluids, and the equation for the non-conserved (granular) 
temperature field is a balance equation for the energy. 
A collisional model for a 2D system was also developed 
on the same basis (see Sec. IV of Ref. (l^); however, 
in this model the stationary temperature is density inde¬ 
pendent and the pressure increases monotonically with 
density, and accordingly there is no two-phase coexis¬ 
tence, contrary to our experimental observation. In this 
model, the granular temperature in a homogeneous sta¬ 
tionary state was assumed to be determined by a balance 
between energy dissipation and injection |45j . but the en¬ 
ergy dissipation due to an effective friction proportional 
to a particle velocity was not considered. Although it is 
interesting to take this effect into account in this theory 
and to study whether it can induce two-phase coexis¬ 
tence, here we take a more phenomenological approach. 

To take this velocity-dependent dissipative force into 
account in a natural manner, we consider a thermal con¬ 
dition to maintain a steady state in a quasi-2D system, 
on the basis of a Langevin-like equation of motion (see, 
e.g., gS]). 

D. A physical principle behind two-phase 
coexistence in a non-equilibrinm steady state. 

Now we consider a principle behind two-phase coex¬ 
istence in a non-equilibrium steady state. For a ther¬ 
modynamic system, in order for two phases of the same 
substance in contact to be in equilibrium, there must 
be mechanical, thermal, and chemical equilibrium, i.e., 
(i) equal pressure, (ii) equal temperature, and (hi) equal 
chemical potential. For our driven granular system, how¬ 
ever, condition (ii) is not be applicable because of the out- 
of-equilibrium nature. On the basis of a physical picture 
we described in Sec. IV.B, here we focus our attention 
on this condition, using a model of granular Brownian 
motion [46] . 

The change in particle velocity due to a binary in¬ 
stantaneous collision between particle i and j with the 
same mass m is given by = v' — [(v' — v') • n] n, 

Vj = v' -I- [(v' — v') • ri] h, where v and v' are the 


velocity after and before the collision respectively, n is 
the unit vector jointing the centres of particles and a is 
the restitution coefficient (0 < a < I) which is equal to 
1 in the elastic case. In order to maintain a steady state, 
we need an external energy source that is coupled to ev¬ 
ery particle in the form of a thermal bath. In our case, 
this external energy is supplied by collisions with vibrat¬ 
ing walls confining a granular monolayer. The motion of 
a particle i is then described by the following stochastic 
equation [46] : 

= fi(t)-7Vi(t)+Cb( 0 - ( 1 ) 

Here ii is the force taking into account the collisions 
with other particles, 7 = m/r is a drag coefficient char¬ 
acterizing the velocity decay towards a steady state, 
whose timescale we express by r, and Cb(^) is random 
white force noise exerted by particle-wall collisions, with 
(Ch(0) = 0 and {Cb,ia{t)Cb,jp{t')) = 2Tb^SijSo.i3S{t - t') 
{Th being the effective bath temperature). Here we note 
that the viscous term in Eq. Q takes into account the 
friction among particles and energy transfers between 
various degrees of freedom (e.g, the friction between par¬ 
ticles and the walls) [51] . 

Here we mention the random nature of force noise. In 
our system, the excitation itself has a well-defined fre¬ 
quency, but the randomization of the energy injection by 
interparticle collision and the wall roughness leads to a 
self generated effective white bath in a long timescale. 
The Langevin-type approach can be rationalized by the 
fact that Brownian dynamics simulations can fully repro¬ 
duce the behavior of a driven monolayer granular system 
of weakly polydisperse steel balls even without any ad¬ 
justable parameters [22]. Its validity to a rubber ball 
system is not obvious, but our visual inspection of the 
motion of steel and rubber balls at least indicates the 
same type of random motion. Although the exact nature 
should be investigated, e.g., by measuring the velocity 
distribution function in the future, this visual inspection 
and the consideration in Sec. VI.C tell us that in a steady 
state the physics may be the same between the two sys¬ 
tems at least on a phenomenological level. So we assume 
that the randomization also takes place efficiently for a 
rubber ball system. 

A stationary state is maintained since the effect of the 
external energy source balances the energy lost by inter¬ 
particle collisions and particle-wall collisions. The key 
parameters of the system is the characteristic velocity 
decay time, t , and the packing fraction (j). Multiplying 
the above stochastic equation Eq. Q by v and averaging 
yields 

= (v(i) • f( 0 ) - 7 (v(t)^) + (v(t) • Cb(i))( 2 ) 

Here we note that the granular temperature T* is de¬ 
fined as T* = and thus the left-hand side of the 

above equation can be written as 2dT*Jdt. The first 
term in the right-hand side of Eq. ([^ can be writ¬ 
ten as (v(t) • f(t)) = -{AE)coh where AE = gH(l - 




12 


[(vi — V 2 ) • fi]^, represents the average energy dissi¬ 
pation rate due to interparticle collisions. (• • • )coi is the 
collision average, whose expression is known for granular 
gas [35] but not for high-density liquid or solid because 
of the difficulty associated with not only excluded vol¬ 
ume effects but also recollisions and memory effects. We 
express this term (v(t) • f(t)) = —{AE)coi as 21). Al¬ 
though its exact form is not clear, it should be a decreas¬ 
ing function of a and an increasing function of (j>, and T*. 
The second term is rewritten from the definition of T* as 
2T*/r. Finally, the third term is the total energy input 
rate and written as (v(t) • Cb{t)) = SyTb/m = 2TbjT. 

Here we briefly consider the role of interparticle fric¬ 
tion. Although the friction coefficient itself may be sim¬ 
ilar between the steel and rubber balls, the rubber ball 
may have a larger coefficient of rolling friction with the 
walls as well as another ball upon collision than the steel 
ball does. On a phenomenological level, this effect can 
be included in the coefficient 7 ; then, 7 , or r, should be 
an increasing function of </). The friction also affects the 
term V (see, e.g.. Ref. |3Z])- We also note that the pres¬ 
ence of interparticle friction also induces forces tangential 
to the domain interface. This may perturb the interface 
position and induce large fluctuations of the interface. 
Such a signature can indeed be seen in Fig. although 
it might be largely due to the large pressure fluctuations 
unique to a driven granular system (see Sec. IV.B). For 
simplicity, however, we do not consider this dynamical 
effect when we discuss the phase coexistence. 

Although the difference in the friction property be¬ 
tween steel and rubber balls may play a role in the ob¬ 
served difference in the phase transition behavior, the 
coefficient of rolling friction is generally small for smooth 
spherical particles. Thus, we argue that the difference in 
the restitution coefficient between the steel and rubber 
ball should be the major source of the difference in the 
energy dissipation between the two systems. 

By incorporating all the above factors, Eq. ([^ can be 
expressed as follow: 


constraint £{(/), T*) = 0. Here we assume that £ is a 
function of the two state variables, (j) and T*. First the 
condition of equal pressure P{(j)L,T^) = P{(j)s,Tg) = Pq 
immediately tells us that r£ > Tg since (pL < i'S: which 
is seen in Fig. id). We note that this condition may be 
robust since it is of mechanical origin, as described before. 
Then the two conditions £{(j)L,T^) = £{(j)s,Tg) = 0, 
together with Eq. tell us that the denser solid 

phase dissipates more energy. Another condition is the 
chemical equilibrium condition, which is also a subtle 
issue. Chemical-balance conditions of two coexisting 
phases with different temperatures should be obtained as 
a steady state solution of the relevant kinetic equations 
also considering bond orientational and translational or¬ 
dering. It may be a promising way to introduce a non¬ 
equilibrium free energy [48] . J^(((), T*, P), and consider 
the balance of dP/dp between the two phases, which are 
expected to be the necessary conditions to maintain a 
steady state. Our discussion is purely phenomenological, 
and furthermore we also need to consider Thus, a 
more rigorous approach is highly desirable in the future. 

The lever rule As/A^ = {p — Pl)/{Ps — P) is then 
obtained from the condition (1) N = Ng + Ns, i.e., 
PlAl + PsAs = pA, where pi, Ni, and Ai, are the par¬ 
ticle area fraction, the number of particles, and the total 
area of phase i, and (2) A = Al + As- The similar lever 
rule was also observed for the monolayer-bilayer transi¬ 
tion (see, e.g., [331ISH])- The crucial difference from a 
thermal system arises from the fact that in our system 
the granular temperature (i.e., the kinetic energy), T*, is 
different between the two phases. The above condition 
£ = 0 is required for maintaining the dynamical steady 
state dissipating energy and it may be this condition that 
plays the most crucial role in the phase selection in a non¬ 
equilibrium state. 


E. Dissipation-induced wetting 


at T T 


(3) 


In a steady state, we should have the relation £ = 0. 
This relation is just a consequence of energy conservation, 
simply implying that the energy input to a 2D granular 
system is partially dissipated by inelastic interparticle 
collisions and by viscous damping due to the particle- 
wall interactions. For an elastic system where a = 1 and 
thus I? = 0 , this relation together with £ = 0 reduces 
to the relation = T*. Strictly speaking, we should 
also consider the energy flux arising from the spatial in¬ 
homogeneity of p and T* (see below), but we tentatively 
ignore it here since it is not relevant to the description 
od a steady state. 

On the basis of the above physical picture, we argue 
that the strong discontinuous nature of the transition of 
an inelastic system found here is a consequence of the in¬ 
homogenization of a system under the above-mentioned 


The side wall-particle collisions dissipate energy with 
a rate different from bulk, which affects the wettability 
of a phase to a solid side wall, as shown in Fig. From 
Eq. we can infer that a wall harder than particles 
tends to wet a liquid phase of high T* rather than a solid 
phase of low T* whereas a wall softer than particles tends 
to wet a solid phase rather than a liquid phase. This is 
simply because particles near a hard/soft side wall tend 
to have higher/lower T*. This spatial inhomogeneity of 
T* near a side wall should be coupled to the local area 
fraction p to satisfy the following steady state condition 
for the balance of the heat flux q |35j : 

q = -kVT* - XVp = 0, (4) 

where k and A are transport coefficients. This condition 
qualitatively explains why the solid phase with a higher 
p tends to wet to the softer wall, which locally lowers T*. 
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F. Dissipation-induced demixing 


Demixing shown in Fig. 11 can be explained by the 
fact that the energy loss upon collision with the top and 
bottom plates is smaller for steel balls than for rubber 
balls. The resulting higher kinetic energy of steel balls 
is a reason why steel balls tends to be located in a liq¬ 
uid phase of hight T*. This is again consistent with the 
physical principle discussed above. A higher kinetic en¬ 
ergy of a steel ball should lead to a larger specific area 
per particle for it. Such a tendency can be seen in Fig. 
[Tl|a) and (b). This can also be explained by the con¬ 
dition of the balance of the heat flux q (see Eq. 0): 
which tells us the negative correlation between the local 
effective temperature T* and the local area fraction (j). 


V. CONCLUSION AND OUTLOOK 


icantly [SH [SSI [SO]- In active matter, and particularly 
in living systems, interactions between active objects are 
often dissipative or non-conservative. Thus, it is quite 
interesting and important to study effects of dissipative 
interactions such as inelasticity and local (or nonlocal) 
friction on the state selection of active matter in a sys¬ 
tematic manner. We hope that our study can aid the 
understanding of fundamental roles of dissipative inter¬ 
actions in self-organization of out-of-equilibrium systems. 
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To summarize, we find that energy dissipation plays a 
crucial role in self-organization of driven monolayer gran¬ 
ular matter, which allows the coexistence of states with 
different effective temperatures, contrary to the phase 
coexistence in a thermal system. The two-phase coexis¬ 
tence obeys the lever rule as in a thermodynamic first- 
order transition; however, the underlying selection rule is 
fundamentally different in the sense that the energy dis¬ 
sipation rate, which is an intrinsically non-equilibrium 
quantity, is the key factor of the phase selection. Since 
our discussion is phenomenological, however, it is desir¬ 
able to theoretically describe the coexistence conditions 
in a non-equilibrium steady state in a more rigorous man¬ 
ner. There are also many fundamental open questions 
such as what determines the interface profile and what is 
the nature of fluctuations of the interface. We also show 
that it is possible to separate particle species by solid¬ 
ification using the inelasticity contrast, which is similar 
to purification of materials including impurities by crys¬ 
tallization in thermodynamic systems. Our findings may 
shed light on a general principle governing the state se¬ 
lection of granular matter far from equilibrium, which 
should also be important for our understanding of indus¬ 
trial processing of granular materials by vibration and 
flow. 

Finally we expect that a similar principle may hold 
for the state selection in active matter, which is another 
important class of out-of-equilibrium systems. For ac¬ 
tive systems, there have recently been many studies on 
a liquid-solid transition glass transition [51H5S] . 

demixing of self-propelled particles [57MD] and rotors 
(saiiiiiig. These studies have elucidated unique char¬ 
acters of the phenomena distinct from their thermody¬ 
namic counterparts. It has been clarified that the cou¬ 
pling between local energy input (or motility) and den¬ 
sity [63] plays a crucial role in the state selection (see, 
e.g., [64||65])- Inclusion of inelastic interactions, nonlo¬ 
cal viscous dissipation, and local friction may also alter 
the nature of the state selection in such a system signif- 


APPENDIX A: COMPARISON OF OUR STUDY 
WITH PREVIOUS WORKS IN QUASI-2D 
SYSTEMS 

A-1: Bilayer-forming solidification in quasi-2D 
driven granular systems 

Here we compare the results of our work with those of 
previous works done in a quasi-2D situation, for which 
bilayer formation is allowed |23H2S]- All the previ¬ 
ous works, which reported liquid-solid coexistence, were 
made for a cell height h of about 1.7-2.0 times of the 
particle diameter d. Accordingly, these systems can form 
bilayers, and the liquid-solid transition in these works al¬ 
ways accompanies monolayer (liquid)-bilayer (solid) tran¬ 
sition. Thus, the situation is essentially different from 
ours, where the cell thickness is thin enough to avoid bi¬ 
layer formation. Interestingly, this extra degree of free¬ 
dom, bilayer formation, makes the physics essentially dif¬ 
ferent from our case in which granular particles always 
form a monolayer. 

The crucial difference can be seen most clearly in the 
dependence of the phase behavior on the strength of in¬ 
elasticity. Urbach and his coworkers reported that inelas¬ 
ticity significantly expands the low-density liquid region 
for a system in which bilayer formation is allowed [inii27]. 
Furthermore, computer simulations showed that the or¬ 
dered phase is not present at any vibration amplitude 
when the inelasticity is large. Similar behaviors were 
also reported by Clerc et al. [H]. These results clearly 
indicate that ordering is suppressed hy inelasticity. We 
stress that this is completely the opposite to our case: In 
our case inelasticity changes the nature of the transition 
from continuous to discontinuous and the upper bound 
density of a disordered liquid state is lower for the in¬ 
elastic system than for the elastic one, i.e., inelasticity 
helps the ordering. Another important difference comes 
from the fact that bilayer formation leads to a change in 
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the number density of the bottom layer because of the 
conservation of the total number of particles. This extra 
degree of freedom also causes a difference in the physics 
between the two types of systems. 

Lobkovsky et al. |27j showed by numerical simu¬ 
lations that the monolayer-bilayer transition becomes 
rather continuous for random forcing. This is an inter¬ 
esting observation in the sense that the way of driving 
affects the nature of the transition. However, they also 
showed that inelasticity suppresses the onset of the or¬ 
dered phase with random forcing, as is observed in the 
vibrating system. We stress that this tendency is again 
the opposite to our case. 

Furthermore, we note that the monolayer-bilayer for¬ 
mation can take place as a function of T [24l [43] , whereas 
our transition cannot be induced by changing T, as we 
can see from the phase diagrams in Fig. A change from 
a discontinuous to a continuous nature of the F-induced 
transition was observed with an increase in the normal¬ 
ized cell thickness h/d from 1.83 to 1.94 [33]. Here the 
key control parameter is the cell thickness, which affects 
the symmetry of the bilayer solid phase (square or hex- 
atic symmetry) [24l[43], and not the degree of inelasticity. 
This fact and the above-mentioned crucial difference in 
the F-dependence of the transition also tell us that the 
physics is essentially different between the two cases. 

The phase ordering behavior in granular matter was 
proposed to be generally expressed by the Ginzburg- 
Landau-type free energy (or potential) and relevant equa¬ 
tions of motion [48]. On the basis of this picture, a theo¬ 
retical explanation for bilayer formation was proposed by 
Clerc et al. [33] . The thermodynamic force on the density 
field u mainly comes from the effective (non-equilibrium) 
pressure gradient calculated from this potential. The 
authors also included the friction term for the u field 
and random force noises. The kinetic equation governing 
the phenomena was shown to be non-diffusive and in¬ 
clude not only du/dt (i.e, dcj)/dt) but also d'^ujdt} (i.e., 
d'^(j)ldt^) [33]. The sum of these forces leads to accelera¬ 
tion of the field u. The reaction-diffusion-type equation 
of motion with these forces captures the travelling wave 
features. This is markedly different from our monolayer 
case, where we have never observed such travelling waves. 

We confirmed by comparing results of experiments and 
those of Brownian dynamics simulations that our driven 
monolayer granular system obeys Langevin (or Brown¬ 
ian) dynamics with viscous damping and random noise 
|22j . This indicates that there is no acceleration term, 
which is proportional to d^(j)/dt^, in our system. Thus, 
the above theory cannot be applied to our monolayer 
case. This presence or absence of the acceleration term 
crucially affects the dynamics and make our monolayer 
system distinct from a bilayer-forming system. Never¬ 
theless, it does not influence a steady state, since there 
dXjdt = 0 for any quantity X, In a steady state, thus, 
the two-phase coexistence is determined by the Ginzburg- 
Landau-type potential [3H] . This coexistence is then con¬ 
trolled solely by the parameter e, which is a coefficient 


of the quadratic term (u^) in the potential. This param¬ 
eter e was assumed to be proportional to the inverse of 
the compressibility coefficient, as in the case of a usual 
thermodynamic gas-liquid coexistence. This (negative) 
compressibility tells us how easily bilayers can be formed. 
In this model, the coexistence is determined by the equal 
pressure and the equal non-equilibrium chemical poten¬ 
tial under the mass conservation. The parameter e was 
shown to be given by the derivative of the momentum 
flux with respect to the density and its negative value, 
which leads to the coexistence, reflects the fact that the 
granular temperature is lower for a higher density [48] . 
Thus, this non-equilibrium chemical potential allows the 
granular temperature to be different between the two co¬ 
existing phases. 

Here we propose an intuitive explanation for the role 
of inelasticity in the liquid-solid transition accompany¬ 
ing bilayer formation, i.e., a relation between the above 
e and inelasticity. For a quasi-2D system which can form 
a bilayer, even for elastic particles there is a liquid-solid 
coexistence (see, e.g.. Fig. 3 in [33]). This is because the 
bilayer formation should be easier for more elastic parti¬ 
cles since a smaller loss of the kinetic energy associated 
with interparticle collisions allows particles in the bot¬ 
tom layer to more efficiently jump to the top layer. With 
an increase in the inelasticity, the discontinuous nature 
becomes weaker and eventually the transition becomes 
continuous for particles having a small enough restitution 
coefficient. This termination of the discontinuous first- 
order-like transition was called a critical point (e = 0) 
by Glerc et al. [33]. Criticality was also observed ex¬ 
perimentally, using F as a control parameter, when F 
approaches a critical value [43] . The theory may be valid 
for a quasi-2D driven granular system where bilayer for¬ 
mation is allowed. However, we emphasize again that the 
above-mentioned dependence on the inelasticity is com¬ 
pletely opposite to what we find in our system: Our cen¬ 
tral finding is that an increase in inelasticity changes the 
nature of the 2D liquid-solid transition from (thermal- 
like) two-step continuous transitions to a one-step dis¬ 
continuous transition. For bilayer formation the kinetic 
energy has to overcome both gravity and energy loss due 
to inelastic collisions. This tells us that more elastic par¬ 
ticles can form a bilayer more easily, as described above. 
To our knowledge, our work is the first to show a discon¬ 
tinuous first-order-like liquid-solid transition of a mono- 
layer granular system in a high F regime, whose nature is 
primarily controlled by energy dissipation, as discussed 
in Sec. IV. 


A-2: Pattern formation due to inhomogeneization of 
the energy input 

Next, we mention another ‘apparently’ similar behav¬ 
ior, but which also has a very different physical origin. 
Olafsen and Urbach [33] reported clustering or ordering 
upon decreasing the vertical vibration amplitude for a 
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quasi-2D system, but which has no upper plate and thus 
is not confined by two parallel walls unlike our case: At 
large F, particle correlations exhibit only short-range or¬ 
der as in the case of equilibrium 2D hard-sphere gases, 
but lowering F cools the system, resulting in a dramatic 
increase in correlations leading to either clustering or an 
ordered state. Further cooling forms a collapse: a con¬ 
densate of motionless balls coexisting with a less dense 
gas. Measured velocity distributions are non-Gaussian, 
showing nearly exponential tails. In our systems, we ob¬ 
served similar phenomena below a critical value of F^, 
which was shown in Fig. |^b) for a rubber system as the 
“inhomogeneous excitation” state (the gray region). For 
a steel ball system this phenomenon was observed for a 
much lower value of F, as described in the main text. We 
note that these phenomena have an essentially different 
physical origin, which is for example discussed in detail 
on the basis of molecular dynamics (MD) simulation [36] 
and a Navier-Stokes granular hydrodynamics m- Nie et 
al. showed by MD simulation |3Bj that at high F the par¬ 
ticle motion is isotropic, and the velocity distributions are 
Gaussian. The deviations from a Gaussian distribution 
at low F is related to the degree of anisotropy in the mo¬ 
tion. Below Fc, the vertical velocity distribution becomes 
bimodal: The cluster particles move with the plate, while 
the gas particles are non-interacting, as they collide pri¬ 
marily with the plate. They proposed that dissipative 
contact forces are responsible for this phenomenon. It 
was also mentioned by Khain and Aranson m that the 
phenomenon can be viewed as a consequence of a nega¬ 
tive compressibility of granular gas. This explanation is 
somewhat similar to that for bilayer formation discussed 
above [29], but the range of F is very different between 
the two: the clustering is observed for low F, whereas 
the ordering accompanying bilayer formation for much 
higher F. In relation to this, it was stated [37] that the 
behavior does not significantly depend on the inelastic¬ 
ity of collisions between the particles; one does not need 
inelastic particle collisions to reproduce experimental ob¬ 
servations. They proposed that the mechanism of phase 
separation occurring at low F is related to the non-trivial 
interplay between the energy injection and the vertical 
temperature of the particles. In any case, thus, an in¬ 
homogeneous granular temperature is primarily a con¬ 
sequence of inhomogeneous energy injection. It was also 
shown that the behavior does not significantly depend on 
the inelasticity of collisions between the particles. This 
clearly indicates that this phenomenon has a physical ori¬ 
gin essentially different from ours. 


A-3: Gas-liquid coexistence in driven inelastic 
particles 

Finally, we mention gas-liquid coexistence observed in 
driven inelastic particles [44]. In this case, the difference 
in the type of particle motion between the two phases was 
found to play a key role in the coexistence: In the dense 


liquid phase, the injected energy is quickly dissipated 
within the bulk by frequent interparticle collisions due 
to a high density. In the dilute gas phase, on the other 
hand, the motion of particles is synchronized with the 
driving, which reduces the relative velocity between par¬ 
ticles and thus the rate of interparticle collisions. Thus, 
two phases with a large difference in the granular tem¬ 
perature coexist. However, this phenomenon is also es¬ 
sentially different from ours, reflecting the difference in 
the cell thickness (many layers vs. monolayer) and the 
nature of the transition (gas-liquid vs. liquid-solid transi¬ 
tion); for example, there is no such synchronized motion 
in our system. 


APPENDIX B: ON THE SPATIAL UNIFORMITY 
OF THE ENERGY INPUT 

Recently Brito et al. [IS] showed that for a granular 
monolayer vibrated between the walls, uniform energy 
input to a system is generally a good assumption. In a 
quasi-2D geometry, in which only a monolayer can exist, 
the system is known to remain homogeneous in the hor¬ 
izontal directions for a wide range of parameters. This 
is due to the presence of a distributed energy source. 
They also pointed out that in the absence of friction, 
this energy source is Galilean invariant and conserves mo¬ 
mentum locally. In such a quasi-2D system, the vertical 
energy scale of grains is fixed by vibration parameters. 
We stress that the energy injection occurs only through 
direct collisions of a particle with the walls, and there 
is no other channel. This means that the energy in¬ 
jection rate is controlled by fpw So, as long as there 
are no direct geometrical restrictions to the vertical mo¬ 
tion of the particles, it is reasonable to assume that the 
energy input is rather homogeneous spatially. We note 
that, according to our observation, there is no overlap of 
particle images projected onto the horizontal plane (see, 
e.g., images in Fig. [^. For steel balls, the input energy 
is homogeneous up to (j) = 0.80, which is much higher 
than the upper bound area fraction of the liquid phase, 
(j)s (=0.695), of the rubber ball system. Importantly, as 
shown in Fig. [^ all the phase transition area fractions, 
(f)h, 4>s, 4 >l, and ())s, are independent of F. This suggests 
that the particle-wall collision frequency fpw may be al¬ 
most the same between the two phases, or rather homo¬ 
geneous spatially, for this range of F. We note that if the 
particle-wall collision frequency fpw (i.e., a vibrational 
parameter) strongly depends on (p, the phase boundary 
compositions should also depend on F. Thus, the energy 
input rate may be assumed to be homogeneous. How¬ 
ever, since the discussion above is qualitative, the spatial 
distribution of the energy input in two-phase coexistence 
needs to be checked carefully by numerical simulations in 
the future. 

For a quasi-2D system, in which a bilayer can be 
formed, the situation is very different. For example, it 
was clearly shown m by mapping the local average rate 
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of energy input that a square phase consisting of bilayer 
corresponds to a region of dramatically reduced energy 
input. However, this is natural since the bilayer forma¬ 
tion itself inevitably accompanies a change in vibrational 


parameters; for example, the effective mass of the vi¬ 
brated object is roughly doubled by bilayer formation. 
We note that such effects are absent for a monolayer sys¬ 
tem, as mentioned above. 
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